Constrained optimization approach to compander design for ofdm papr reduction

ABSTRACT

In the method for constrained optimization for compander design for OFDM PAPR, wherein the improvement comprises the step of converting a Rayleigh amplitude distribution from which superior performing companders are derived, whereby a constrained optimization problem may be solved using Lagrange multipliers.

RELATED APPLICATIONS

This application claims the benefit of U.S. Provisional Application No.62/036,388, filed Aug. 12, 2014. This application is herein incorporatedby reference in its entirety for all purposes.

FIELD OF THE INVENTION

The invention relates to signal processing and, more particularly, tocompanders for reducing peak to average power ratio (PAPR) in orthogonalfrequency division multiplexing (OFDM) signals.

BACKGROUND OF THE INVENTION

Orthogonal frequency-division multiplexing (OFDM) is a method ofencoding digital data on multiple carrier frequencies. OFDM hasdeveloped into a popular scheme for wideband digital communication, usedin applications such as digital television and audio broadcasting, DSLInternet access, wireless networks, powerline networks, and 4G mobilecommunications.

For a sufficiently large number of tones, OFDM signal values areGaussian, or normally, distributed, by the central limit theorem, whichstates that, given certain conditions, the arithmetic mean of asufficiently large number of iterates of independent random variables,each with a well-defined expected value and well-defined variance, willbe approximately normally distributed, regardless of the underlyingdistribution.

The Gaussian distributed signal components generate Rayleigh distributedsignal amplitude values, Rayleigh distributions being continuousprobability distributions for positive-valued random variables,specifically the distribution of the magnitude of a two-dimensionalrandom vector whose coordinates are independent, identicallydistributed, mean 0 normal variables. In signals typical of the currentstate of the art, the Rayleigh distribution is typically heavy-tailed,i.e. a larger proportion of its population rests within its tail thanwould under a normal distribution. The long tail of the Rayleighdistribution, i.e. the portion of the distribution having a large numberof occurrences far from the “head” or central part of the distribution,indicates the presence of large amplitude values, resulting fromconstructive interference. These large amplitude values can saturate thepower amplifier, used in the processing chain, thus creating signaldistortion which reduces demodulation performance.

This causes existing techniques for Orthogonal Frequency-DivisionMultiplexing (OFDM) to suffer from large peak-to-average power ratio(PAPR). The PAPR is the peak amplitude squared (giving the peak power)divided by the mean square value squared (giving the average power).

Many techniques have been developed to mitigate the PAPR problem,including signal companding. Companders are attractive because of theirsimplicity and effectiveness. Companders apply an amplitude weight tothe signal, generally downweighting large amplitude values, whileupweighting smaller values to maintain constant power. The largeamplitude downweighting, while simultaneously maintaining the same powerlevel, reduces the PAPR. The amplitude weighting introduces signaldistortion, so it is important to design companders that introduceminimal distortion to limit demodulation errors, while also providingout-of-band power rejection.

Early solutions adapted companders from speech and audio coding, such asμ-law and A-law companders. Subsequent approaches decomposed theamplitude interval into disjoint regions over which different amplitudeweighting functions were used. Modern companders are derived bytransforming the Rayleigh amplitude distribution into a more favorabledistribution to reduce PAPR. Among these transformations, piecewisetransformations, i.e. a function composed of straight-line sections,have been used, either uniform, linear, nonlinear, or piecewise linear.

A need still exists, however, for solutions which minimally alter theoriginal Rayleigh amplitude distribution.

SUMMARY OF THE INVENTION

One embodiment of the present invention provides a signal transmittingapparatus for an orthogonal frequency division multiplexing (OFDM)system, comprising: a compander configured to compress and expand atransmitted signal wherein the compander comprises a Rayleigh amplitudedistribution to a constrained optimization problem which may be solvedusing Lagrange multipliers.

One embodiment of the present invention provides a signal transmittingapparatus for orthogonal frequency division multiplexing (OFDM) system,comprising: a companding feedback module configured to compress andexpand a transmitted signal; wherein the companding feedback modulecomprises a Rayleigh amplitude distribution to a constrainedoptimization problem which may be solved using Lagrange multipliers.

One embodiment of the present invention provides a method of deriving acompander for orthogonal frequency division multiplexing comprising:converting a Rayleigh amplitude distribution to a constrainedoptimization problem which may be solved using Lagrange multipliers.

Another embodiment of the present invention provides such a methodwherein the Rayleigh amplitude distribution is decomposed into disjointregions over which a plurality of amplitude weighting functions areused.

A further embodiment of the present invention provides such a methodwherein the constrained optimization approach comprises optimizing anobjective function with respect to changeable elements in the presenceof constraints on those elements.

Yet another embodiment of the present invention provides such a methodwherein the Rayleigh amplitude distribution is a probabilitydistribution function and a first the constraint is that the Rayleighamplitude distribution must integrate to unity.

A yet further embodiment of the present invention provides such a methodwherein a second the constraint is that power must be constant acrossthe distribution.

Still another embodiment of the present invention provides such a methodwherein an optimal solution is minimally perturbed to generate onlynon-negative solutions.

A still further embodiment of the present invention provides such amethod of claim 5 wherein the Rayleigh amplitude distribution is definedas:

${f_{R}(x)} = {\frac{2x}{\sigma^{2}}{^{- \frac{x^{2}}{\sigma^{2\;}}}.}}$

Even another embodiment of the present invention provides such a methodwherein further comprising determining a function g(x) which minimizes:

∫₀^(A)(f_(R)(x) − g(x))²x.

An even further embodiment of the present invention provides such amethod wherein g(x) is chosen to be of piecewise linear form.

A still even another embodiment of the present invention provides such amethod wherein g(x) is defined as

${g(x)} = {\sum\limits_{i = 1}^{N}{{U_{i}(x)}.}}$

A still even further embodiment of the present invention provides such amethod wherein the dynamics of the system are defined as: Λ(β₁, . . .β_(N+1), λ₁, λ₂)=F(β₁, . . . β_(N+1))+λ₁[g₁(β₁, . . .β_(N+1))−1]+λ₂[g₂(β₁, . . . β_(N+1))−σ²].

Still yet another embodiment of the present invention provides such amethod wherein a first constraint, that the cumulative distributionfunction must integrate to unity, is expressed as: g₁(β₁, . . .β_(N+1))=1.

A still yet further embodiment of the present invention provides such amethod wherein a second constraint, that power must be constant, isexpressed as: g₂(β₁, . . . β_(N+1))=σ².

Even yet another embodiment of the present invention provides such amethod wherein the constrained optimization problem can be defined as:Minimize:

F(β₁, …  β_(N + 1)) = ∫₀^(A)(f_(R)(x) − g(x))²x

for β_(i), subject to

∫₀^(A)g(x)x = 1  and  ∫₀^(A)x²g(x)x = ∫₀^(∞)x²f_(R)(x)x = σ²,

using Lagrange multipliers.

An even yet further embodiment of the present invention provides such amethod wherein the compander is defined as:

${F_{y}^{- 1}{F_{R}(x)}} = {\sqrt[ \pm ]{\frac{1 - ^{- \frac{x^{2}}{\sigma^{2}}} - Z_{i - 1}}{\gamma} - \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)} - \frac{\delta}{2\gamma}}$

for a_(i)≦x<a_(i+1), wherein

$Z_{i} = {Z_{i - 1} + {\frac{{\beta_{i + 1}\left( {a_{i + 1} - a_{i}} \right)} + {\beta_{i}\left( {a_{i + 1} - a_{i}} \right)}}{2}.}}$

Still even yet another embodiment of the present invention provides sucha method wherein a decompander for use with the compander is defined as:

${x = {{{sgn}(x)}\sqrt{{- \sigma^{2}}{\ln \left( {{- {\gamma \left( {y + \frac{\delta}{2\gamma}} \right)}^{2}} - e + \frac{\delta^{2}}{4\gamma} + 1 - Z_{i - 1}} \right)}}}},$

where the sgn function refers to retaining either the positive ornegative square-root, depending on the sign of x; and wherein the rangeof the compander, when a_(i)≦x<a_(i+1), then y_(MIN)≦y<y_(MAX) where:y_(MIN)=min{y(a_(i)),y(a_(i+1))}, y_(MAX)=max{y(a_(i)),y(a_(i+1))}, canbe defined as

${y\left( a_{i + 1} \right)} = {\sqrt[ \pm ]{\frac{1 - ^{- \frac{a_{i + 1}^{2}}{\sigma^{2\;}}} - Z_{i - 1}}{\gamma} - \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)} - {\frac{\delta}{2\gamma}.}}$

A still even yet further embodiment of the present invention providessuch a method wherein a Hessian Matrix is used to determine criticalpoint types.

The features and advantages described herein are not all-inclusive and,in particular, many additional features and advantages will be apparentto one of ordinary skill in the art in view of the drawings,specification, and claims. Moreover, it should be noted that thelanguage used in the specification has been principally selected forreadability and instructional purposes, and not to limit the scope ofthe inventive subject matter.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a graph showing piecewise linear approximation to a preferredembodiment of the invention;

FIG. 2 is a graph showing optimal solution, A=1.1;

FIG. 3 is a graph showing optimal solution, A=1.2;

FIG. 4 is a graph showing optimal solution, A=1.3;

FIG. 5 is a graph showing optimal solution, A=1.4;

FIG. 6 is a graph showing optimal solution, A=1.5;

FIG. 7 is a graph showing optimal solution, A=1.6;

FIG. 8 is a graph showing optimal solution, A=1.7;

FIG. 9 is a graph showing optimal solution, A=1.8;

FIG. 10 is a graph showing positive perturbed solution; A=1.4;

FIG. 11 is a graph showing positive perturbed function;

FIG. 12 is a graph showing perturbation solution; A=1.4, k=5;

FIG. 13 is a graph showing perturbation solution; A=1.4, k=6;

FIG. 14 is a graph showing perturbation solution; A=1.4, k=7;

FIG. 15 is a graph showing perturbation solution; A=1.4, k=8;

FIG. 16 is a graph showing symbol error rate performance;

FIG. 17 is a graph showing PAPR reduction performance;

FIG. 18 is a graph showing out-of-band power rejection; and

FIG. 19 is a graph showing sample Lagrange Companders.

DETAILED DESCRIPTION

The present disclosure describes new companders that use a constrainedoptimization approach for reducing peak-to-average power ratio (PAPR) inorthogonal frequency-division multiplexing (OFDM) signals. A constrainedoptimization approach, in this context, involves optimizing an objectivefunction, specifically an energy function, with respect to changeableelements, i.e. variables, in the presence of constraints on thosevariables.

Companders, devices which use both compression and expansion to improvedynamic range and signal-to-noise ratio, in accordance with the presentdisclosure provide symbol error rate performance improvements overcurrent state-of-the art companders. In digital communications, symbolrate (also known as baud or modulation rate) is the number of symbolchanges (waveform changes or signalling events) made to the transmissionmedium per second using a digitally modulated signal. The symbol rate ismeasured in baud (Bd) or symbols per second. Each symbol can representor convey one or several bits of data. The symbol error rate is simplythe rate at which signalling events fail to convey the intended data.

The newly-designed companders also provide design flexibility, therebyexpanding the space of tradeoffs between demodulation performance, PAPRreduction, and out-of-band power rejection. Furthermore, the newcompanders provide solutions in operating regions where certain currentcompanders fail to exit; solutions may be derived for cutoff amplitudevalues that are unobtainable using other companders. The new compandersolutions enable tuning the cutoff amplitude value based on poweramplifier bandwidth.

Below, we formulate the constrained optimization problem and derive thecompander and decompander forms. Through numerical simulation, wegenerate performance results demonstrating the capability of the newcompanders.

1. COMPANDER DESIGN USING CONSTRAINED OPTIMIZATION

This section describes the constrained optimization problem andsolution, and derives the compander and decompander.

1.1. Constrained Optimization Problem Definition

Denote the Rayleigh amplitude distribution by ƒ_(R)(x) for:

${f_{R}(x)} = {\frac{2x}{\sigma^{2}}^{- \frac{x^{2}}{\sigma^{2}}}}$

In the above equation, ƒ_(R)(x) represents the probability densityfunction of the amplitude values of the OFDM signal. This is astatistical description of the spread of amplitude values of the OFDMsignal, which follow a Rayleigh distribution. We then look for afunction g(x)which minimizes:

$\begin{matrix}{\int_{0}^{A}{\left( {{f_{R}(x)} - {g(x)}} \right)^{2}{x}}} & (1)\end{matrix}$

This function is used because we are making the assumption, or operatingunder the hypothesis, that minimal distortion of the natural signaldistribution (i.e., the Rayleigh distribution) will incur minimal signaldistortion of the companded function and hence reduce the demodulationerrors (i.e., reduce the symbol error rate).

We then choose g(x) to be of piecewise linear form:

${g(x)} = {\sum\limits_{i = 1}^{N}{U_{i}(x)}}$

In the above equation, the piecewise linear form was chosen since itextends previous work on using single linear components, with the viewthat the increased flexibility of the piecewise approach will lead tobetter performing companders. Another reason it was used is that itallows a closed form expression for the compander weights to be found,by solving a linear system of equations.

Still referring to the above equation, each U_(i)(x) is acompactly-supported linear segment, in this context, meaning a linesegment of finite extent, i.e., of finite length. The U_(i)(x) isdefined by equation (2) and for each index i, represents the equationfor one of the linear segments shown in FIG. 1.

$\begin{matrix}{{U_{i}(x)} = \left\{ \begin{matrix}{{\left( \frac{\beta_{i + 1} - \beta_{i}}{a_{i + 1} - a_{i}} \right)\left( {x - a_{i\;}} \right)} + \beta_{i}} & {{{for}\mspace{14mu} x} \in \left\lbrack {a_{i},a_{i + 1}} \right\rbrack} \\0 & {otherwise}\end{matrix} \right.} & (2)\end{matrix}$

and (a₁, . . . a_(N+1)) is a partition of the interval [0, A] and theβ_(i) are the ordinates at x=a_(i) (FIG. 1). This interval is a standardinterval used in PAPR reduction techniques, and represents an intervalof amplitude values, covering amplitude values between 0 and A. Anyamplitude value in the original OFDM signal that is larger than A getsmapped to something less than or equal to A, so A is the maximum valuethe companded value can have.

With the assumed form of g(x), we want to minimize:

$\begin{matrix}{{F\left( {\beta_{1},{\ldots \mspace{14mu} \beta_{N + 1}}} \right)} = {\int_{0}^{A}{\left( {{f_{R}(x)} - {g(x)}} \right)^{2}\ {x}}}} & (3)\end{matrix}$

as a function of the β_(i). The function that we are minimizingrepresents a measure of the deviation of the new, companded functionamplitude value distribution from the original OFDM Rayleigh amplitudedistribution. The beta values are the values of the piecewise linearcomponent functions at the endpoints, as shown in FIG. 1. (i.e., theordinates), and described by the functions U_(i)(x).

There are two constraints on g(x); g(x) is a probability densityfunction and so must integrate to unity value over xε[0, A]. This isbecause, by definition, a probability density function gives an amountof probability per unit of something, in this case, per unit of signalamplitude value. Since the total probability of all of the possibilitiesadds up (integrates up) to one, then the integral of the probabilitydensity function over all of its values must equal the total probabilityvalue of one.

$\begin{matrix}{{\int_{0}^{A}{{g(x)}\ {x}}} = 1} & (4)\end{matrix}$

This is the unity cumulative distribution function (c.d.f.) constraint.The c.d.f. itself describes the probability that a real-valued randomvariable x with a given probability distribution will be found to have avalue less than or equal to x. In the case of a continuous distribution,it gives the area under the probability density function (p.d.f.) fromminus infinity to x.

The second constraint is the constant power constraint:

$\begin{matrix}{{\int_{0}^{A}{x^{2}{g(x)}\ {x}}} = {{\int_{0}^{\infty}{x^{2}{f_{R}(x)}\ {x}}} = \sigma^{2}}} & (5)\end{matrix}$

In the above equation, equation 5, the integral on the right representsthe total power across all possible amplitude values, and its value isequal to the constant value sigma times sigma, hence the total power isthe constant value sigma times sigma.

There is actually a third constraint on g(x); because g(x) is a p.d.f.,we have a non-negativity constraint:

g(x)≧0  (6)

However, we forgo the non-negativity constraint, because, as isdiscussed later, we can slightly perturb the optimal solution togenerate non-negative solutions. The constrained optimization problemthen becomes: minimize (3) for the β_(i) subject to (4) and (5). Wesolve this using Lagrange multipliers, a strategy for finding the localmaxima and minima of a function subject to equality constraints.

1.2. Calculation of Partial Derivatives

To solve the constrained optimization problem, we need to calculate thepartial derivatives of F. First consider:

$\begin{matrix}{{F\left( {\beta_{1},{\ldots \mspace{14mu} \beta_{N + 1}}} \right)} = {\int_{0}^{A}{\left( {{f_{R}(x)} - {\sum\limits_{i = 1}^{N}{U_{i}(x)}}} \right)^{2}\ {x}}}} & (7)\end{matrix}$

This equation results from eq. (1), where we use the piecewise linearform for function g, and call the whole thing F, and in F we explicitlyidentify the fact that the expression is now a function of the betaparameters, as shown from eq 2.

Because the integrand in (7) is continuously differentiable andintegrable for x≠β_(i), then differentiation is permitted under theintegral sign and:

$\begin{matrix}\begin{matrix}{\frac{\partial F}{\partial\beta_{j}} = {\int_{0}^{A}{\frac{\partial}{\partial\beta_{j}}\left( {{f_{R}(x)} - {\sum\limits_{i = 1}^{N}{U_{i}(x)}}} \right)^{2}\ {x}}}} \\{= {\int_{0}^{A}{2\left( {{f_{R}(x)} - {\sum\limits_{i = 1}^{N}{U_{i}(x)}}} \right)\left( {{- \ \frac{\partial}{\partial\beta_{j}}}{\sum\limits_{i = 1}^{N}{U_{i}(x)}}} \right){x}}}} \\{= {\int_{0}^{A}{2\left( {{f_{R}(x)} - {\sum\limits_{i = 1}^{N}{U_{i}(x)}}} \right)\left( {- \ {\sum\limits_{i = 1}^{N}\frac{\partial{U_{i}(x)}}{\partial\beta_{j}}}} \right){x}}}} \\{= {{2{\int_{0}^{A}{{- {f_{R}(x)}}{\sum\limits_{i = 1}^{N}\frac{\partial{U_{i}(x)}}{\partial\beta_{j}}}}}} + {\sum\limits_{k = 1}^{N}{{U_{k}(x)}{\sum\limits_{i = 1}^{N}{\frac{\partial{U_{i}(x)}}{\partial\beta_{j}}{x}}}}}}} \\{= {{{- 2}{\int_{0}^{A}{{f_{R}(x)}{\sum\limits_{i = 1}^{N}{\frac{\partial{U_{i}(x)}}{\partial\beta_{j}}{x}}}}}} +}} \\{{2{\int_{0}^{A}{\sum\limits_{k = 1}^{N}{{U_{k}(x)}{\sum\limits_{i = 1}^{N}{\frac{\partial{U_{i}(x)}}{\partial\beta_{j}}{x}}}}}}}} \\{= {{{- 2}{\sum\limits_{i = 1}^{N}{\int_{a_{i}}^{a_{i + 1}}{{f_{R}(x)}\frac{\partial{U_{i}(x)}}{\partial\beta_{j}}{x}}}}} +}} \\{{2{\sum\limits_{i = 1}^{N}{\int_{0}^{A}{{U_{i}(x)}\frac{\partial{U_{i}(x)}}{\partial\beta_{j}}{x}}}}}} \\{= {{{- 2}{\sum\limits_{i = 1}^{N}{\int_{a_{i}}^{a_{i + 1}}{{f_{R}(x)}\frac{\partial{U_{i}(x)}}{\partial\beta_{j}}{x}}}}} +}} \\{{2{\sum\limits_{i = 1}^{N}{\int_{a_{i}}^{a_{i + 1}}{{U_{i}(x)}\frac{\partial{U_{i}(x)}}{\partial\beta_{j}}{x}}}}}}\end{matrix} & (8)\end{matrix}$

Now from (2), U_(i)(x) can be rewritten as:

${U_{i}(x)} = {{\left( \frac{x - a_{i}}{a_{i + 1} - a_{i}} \right)\beta_{i + 1}} + {\left( {1 - \frac{x - a_{i}}{a_{i + 1} - a_{i}}} \right)\beta_{i}}}$

so that:

$\begin{matrix}{\frac{\partial U_{i}}{\partial\beta_{j}} = {{\left( \frac{x - a_{i}}{a_{i + 1} - a_{i}} \right)\frac{\partial\beta_{i + 1}}{\partial\beta_{j}}} + {\left( {1 - \frac{x - a_{i}}{a_{i + 1} - a_{i}}} \right)\frac{\partial\beta_{i}}{\partial\beta_{j}}}}} & (9)\end{matrix}$

In the above equation, j represents an index on the beta parameters. Forj=1, equation (9), once simplified, becomes:

$\begin{matrix}{\frac{\partial U_{i}}{\partial\beta_{j}} = \left\{ \begin{matrix}{1 - \frac{x - a_{1}}{a_{2} - a_{1}}} & {{{if}\mspace{14mu} i} = 1} \\0 & {otherwise}\end{matrix} \right.} & (10)\end{matrix}$

For 1<j≦N, (9) becomes

$\begin{matrix}{\frac{\partial U_{i}}{\partial\beta_{j}} = \left\{ \begin{matrix}{1 - \frac{x - a_{j}}{a_{j + 1} - a_{j}}} & {{{if}\mspace{14mu} i} = j} \\\frac{x - a_{j - 1}}{a_{j} - a_{j - 1}} & {{{if}\mspace{14mu} i} = {j - 1}} \\0 & {otherwise}\end{matrix} \right.} & (11)\end{matrix}$

For j=N+1, (9) becomes

$\begin{matrix}{\frac{\partial U_{i}}{\partial\beta_{j}} = \left\{ {\begin{matrix}\frac{x - a_{N}}{a_{N + 1} - a_{N}} & {{{if}\mspace{14mu} i} = N} \\0 & {otherwise}\end{matrix}.} \right.} & (12)\end{matrix}$

So, from (8), and (10), for j=1

$\begin{matrix}{\frac{\partial F}{\partial\beta_{j}} = {{{- 2}{\int_{a_{1}}^{a_{2}}{{f_{R}(x)}\left( {1 - \frac{x - a_{1}}{a_{2} - a_{1}}} \right)\ {x}}}} + {2{\int_{a_{1}}^{a_{2}}{{U_{1}(x)}\ \left( {1 - \frac{x - a_{1}}{a_{2} - a_{1}}} \right){x}}}}}} & (13)\end{matrix}$

From (8), and (11), for 1<j≦N

$\begin{matrix}{\frac{\partial F}{\partial\beta_{j}} = {{{- 2}{\int_{a_{j - 1}}^{a_{j}}{{f_{R}(x)}\left( \frac{x - a_{j - 1}}{a_{j} - a_{j - 1}} \right)\ {x}}}} - {2{\int_{a_{j}}^{a_{j + 1}}{{f_{R}(x)}\ \left( {1 - \frac{x - a_{j}}{a_{j + 1} - a_{j}}} \right){x}}}} + {2{\int_{a_{j - 1}}^{a_{j}}{{U_{j - 1}(x)}\left( \frac{x - a_{j - 1}}{a_{j} - a_{j - 1}} \right)\ {x}}}} + {2{\int_{a_{j}}^{a_{j + 1}}{{U_{j}(x)}\ \left( {1 - \frac{x - a_{j}}{a_{j + 1} - a_{j}}} \right){x}}}}}} & (14)\end{matrix}$

From (8), and (12), for j=N+1

$\begin{matrix}{\frac{\partial F}{\partial\beta_{N + 1}} = {{{- 2}{\int_{a_{N}}^{a_{N + 1}}{{f_{R}(x)}\left( \frac{x - a_{N}}{a_{N + 1} - a_{N}} \right)\ {x}}}} + {2{\int_{a_{N}}^{a_{N + 1}}{{U_{N}(x)}\ \left( \frac{x - a_{N}}{a_{N + 1} - a_{N}} \right){x}}}}}} & (15)\end{matrix}$

Inspection of (13)-(15) show that the following integral forms areneeded:

${K_{R}^{1}\left( {c,d} \right)} = {\int\limits_{c}^{d}{{f_{R}(x)}{x}}}$${K_{R}^{2}\left( {c,d} \right)} = {\int\limits_{c}^{d}{{{xf}_{R}(x)}{x}}}$${K_{U}^{1}\left( {c,d} \right)} = {\int\limits_{c}^{d}{{U_{j}(x)}{x}}}$${K_{U}^{2}\left( {c,d} \right)} = {\int\limits_{c}^{d}{{{xU}_{j}(x)}{x}}}$

In the above equations, K is used as a compact way to represent theintegral on the right side of the equation. After much simplification,straightforward integration produces:

$\begin{matrix}{\mspace{79mu} {{\int\limits_{c}^{d}{{f_{R}(x)}{x}}} = {^{- \frac{c^{2}}{\sigma^{2}}} - ^{- \frac{d^{2}}{\sigma^{2}}}}}} & (16) \\{\mspace{79mu} {{\int\limits_{c}^{d}{{{xf}_{R}(x)}{x}}} = {{c\; ^{- \frac{c^{2}}{\sigma^{2}}}} - {d\; ^{- \frac{d^{2}}{\sigma^{2}}}} + {\frac{\sigma \sqrt{\pi}}{2}\left\lbrack {{\Phi \left( \frac{d}{\sigma} \right)} - {\Phi \left( \frac{c}{\sigma} \right)}} \right\rbrack}}}} & (17) \\{{\int\limits_{c}^{d}{{U_{j}(x)}{x}}} = {{\left\lbrack {\frac{a_{j}\left( {d - c} \right)}{a_{j + 1} - a_{j}} - \frac{\frac{d^{2}}{2} - \frac{c^{2}}{2}}{a_{j + 1} - a_{j}} + \left( {d - c} \right)} \right\rbrack \beta_{j}} + {\left\lbrack {\frac{\frac{d^{2}}{2} - \frac{c^{2}}{2}}{a_{j + 1} - a_{j}} - \frac{a_{j}\left( {d - c} \right)}{a_{j + 1} - a_{j}}} \right\rbrack \beta_{j + 1}}}} & (18) \\{{\int\limits_{c}^{d}{{{xU}_{j}(x)}{x}}} = {{\left\lbrack {\frac{a_{j}\left( {\frac{d^{2}}{2} - \frac{c^{2}}{2}} \right)}{a_{j + 1} - a_{j}} - \frac{\frac{d^{3}}{3} - \frac{c^{3}}{3}}{a_{j + 1} - a_{j}} + \left( {\frac{d^{2}}{2} - \frac{c^{2}}{2}} \right)} \right\rbrack \beta_{j}} + {\left\lbrack {\frac{\frac{d^{3}}{3} - \frac{c^{3}}{3}}{a_{j + 1} - a_{j}} - \frac{a_{j}\left( {\frac{d^{2}}{2} - \frac{c^{2}}{2}} \right)}{a_{j + 1} - a_{j}}} \right\rbrack \beta_{j + 1}}}} & (19)\end{matrix}$

where Φ(x) denotes the Gaussian error function, a special function(non-elementary) of sigmoid shape that occurs in probability,statistics, and other disciplines:

${\Phi (x)} = {\int\limits_{0}^{x}{\frac{2}{\sqrt{\pi}}^{- \xi^{2}}{\xi}}}$

From (18) and (19) we have:

$\begin{matrix}\begin{matrix}{K_{U}^{1} = {{\left\lbrack {\frac{a_{j}\left( {d - c} \right)}{a_{j + 1} - a_{j}} - \frac{\frac{d^{2}}{2} - \frac{c^{2}}{2}}{a_{j + 1} - a_{j}} + \left( {d - c} \right)} \right\rbrack \beta_{j}} +}} \\{{\left\lbrack {\frac{\frac{d^{2}}{2} - \frac{c^{2}}{2}}{a_{j + 1} - a_{j}} - \frac{a_{j}\left( {d - c} \right)}{a_{j + 1} - a_{j}}} \right\rbrack \beta_{j + 1}}} \\{= {{{M_{j}^{U\; 1}\left( {a_{j},a_{j + 1}} \right)}\beta_{j}} + {{N_{j}^{U\; 1}\left( {a_{j},a_{j + 1}} \right)}\beta_{j + 1}}}}\end{matrix} & (20) \\\begin{matrix}{K_{U}^{2} = {{\left\lbrack {\frac{a_{j}\left( {\frac{d^{2}}{2} - \frac{c^{2}}{2}} \right)}{a_{j + 1} - a_{j}} - \frac{\frac{d^{3}}{3} - \frac{c^{3}}{3}}{a_{j + 1} - a_{j}} + \left( {\frac{d^{2}}{2} - \frac{c^{2}}{2}} \right)} \right\rbrack \beta_{j}} +}} \\{{\left\lbrack {\frac{\frac{d^{3}}{3} - \frac{c^{3}}{3}}{a_{j + 1} - a_{j}} - \frac{a_{j}\left( {\frac{d^{2}}{2} - \frac{c^{2}}{2}} \right)}{a_{j + 1} - a_{j}}} \right\rbrack \beta_{j + 1}}} \\{= {{{M_{j}^{U\; 2}\left( {a_{j},a_{j + 1}} \right)}\beta_{j}} + {{N_{j}^{U\; 2}\left( {a_{j},a_{j + 1}} \right)}\beta_{j + 1}}}}\end{matrix} & (21)\end{matrix}$

Using (16), (17), (20), and (21) in (13) through (15) after muchsimplification produces for j=1:

$\begin{matrix}{\frac{\partial F}{\partial\beta_{1}} = {{P_{1}\beta_{1}} + {P_{2}\beta_{2}} + P_{0}}} & (22)\end{matrix}$

Where:

$P_{1} = {{2{M_{1}^{U\; 1}\left( {a_{1},a_{2}} \right)}} - {\left( \frac{2}{a_{2} - a_{1}} \right){M_{1}^{U\; 2}\left( {a_{1},a_{2}} \right)}} + {\left( \frac{2a_{1}}{a_{2} - a_{1}} \right){M_{1}^{U\; 1}\left( {a_{1},a_{2}} \right)}}}$$P_{2} = {{2{N_{1}^{U\; 1}\left( {a_{1},a_{2}} \right)}} - {\left( \frac{2}{a_{2} - a_{1}} \right){N_{1}^{U\; 2}\left( {a_{1},a_{2}} \right)}} + {\left( \frac{2a_{1}}{a_{2} - a_{1}} \right){N_{1}^{U\; 1}\left( {a_{1},a_{2}} \right)}}}$$\mspace{79mu} {P_{0} = {{{- 2}{K_{R}^{1}\left( {a_{1},a_{2}} \right)}} - {\left( \frac{2}{a_{2} - a_{1}} \right){K_{R}^{2}\left( {a_{1},a_{2}} \right)}} - {\left( \frac{2a_{1}}{a_{2} - a_{1}} \right){{K_{R}^{1}\left( {a_{1},a_{2}} \right)}.}}}}$

In equation 22, P is used as a convenient and compact way to representthe expressions shown above.

For 1<j≦N:

$\frac{\partial F}{\partial\beta_{j}} = {{Q_{0}(j)} + {{Q_{1}(j)}\beta_{j - 1}} + {{Q_{2}(j)}\beta_{j}} + {{Q_{3}(j)}\beta_{j + 1}}}$

Where:

${Q_{0}(j)} = {{\left( \frac{- 2}{a_{j} - a_{j - 1}} \right){K_{R}^{2}\left( {a_{j - 1},a_{j}} \right)}} + {\left( \frac{2a_{j - 1}}{a_{j} - a_{j - 1}} \right){K_{R}^{1}\left( {a_{j - 1},a_{j}} \right)}} - {2{K_{R}^{1}\left( {a_{j},a_{j + 1}} \right)}} + {\left( \frac{2}{a_{j + 1} - a_{j}} \right){K_{R}^{2}\left( {a_{j},a_{j + 1}} \right)}} - {\left( \frac{2a_{j}}{a_{j + 1} - a_{j}} \right){K_{R}^{1}\left( {a_{j},a_{j + 1}} \right)}}}$$\mspace{79mu} {{Q_{1}(j)} = {\frac{2{M_{j - 1}^{U\; 2}\left( {a_{j - 1},a_{j}} \right)}}{a_{j} - a_{j - 1}} - \frac{2a_{j - 1}{M_{j - 1}^{U\; 1}\left( {a_{j - 1},a_{j}} \right)}}{a_{j} - a_{j - 1}}}}$${Q_{2}(j)} = {\frac{2{N_{j - 1}^{U\; 2}\left( {a_{j - 1},a_{j}} \right)}}{a_{j} - a_{j - 1}} - \frac{2a_{j - 1}{N_{j - 1}^{U\; 1}\left( {a_{j - 1},a_{j}} \right)}}{a_{j} - a_{j - 1}} + {2{M_{j}^{U\; 1}\left( {a_{j},a_{j + 1}} \right)}} - \frac{2{M_{j}^{U\; 2}\left( {a_{j},a_{j + 1}} \right)}}{a_{j + 1} - a_{j}} + \frac{2a_{j}{M_{j}^{U\; 1}\left( {a_{j},a_{j + 1}} \right)}}{a_{j + 1} - a_{j}}}$$\mspace{79mu} {{Q_{3}(j)} = {{2{N_{j}^{U\; 1}\left( {a_{j},a_{j + 1}} \right)}} - \frac{2{N_{j}^{U\; 2}\left( {a_{j},a_{j + 1}} \right)}}{a_{j + 1} - a_{j}} + \frac{2a_{j}{N_{j}^{U\; 1}\left( {a_{j},a_{j + 1}} \right)}}{a_{j + 1} - a_{j}}}}$

In the equation following paragraph [0071], Q is used as a convenientand compact way to represent the expressions shown above.

M, as used above, is similarly a convenient, compact expression, seen ineq. 20. It is the term inside the first set of brackets in eq. 20.

N, as used above, is also a convenient, compact expression, seen in eq.20. It is the term inside the second set of brackets in eq. 20.

For j=N+1:

$\frac{\partial F}{\partial\beta_{N + 1}} = {{R_{1}\beta_{N}} + {R_{2}\beta_{N + 1}} + R_{0}}$

Where:

$R_{0} = {\frac{{- 2}{K_{R}^{2}\left( {a_{N},a_{N + 1}} \right)}}{a_{N + 1} - a_{N}} + \frac{2a_{N}{K_{R}^{1}\left( {a_{N},a_{N + 1}} \right)}}{a_{N + 1} - a_{N}}}$$R_{1} = {{\left( \frac{2}{a_{N + 1} - a_{N}} \right){M_{N}^{U\; 2}\left( {a_{N},a_{N + 1}} \right)}} - {\left( \frac{2a_{N}}{a_{N + 1} - a_{N}} \right){M_{N}^{U\; 1}\left( {a_{N},a_{N + 1}} \right)}}}$$R_{2} = {{\left( \frac{2}{a_{N + 1} - a_{N}} \right){N_{N}^{U\; 2}\left( {a_{N},a_{N + 1}} \right)}} - {\left( \frac{2a_{N}}{a_{N + 1} - a_{N}} \right){N_{N}^{U\; 1}\left( {a_{N},a_{N + 1}} \right)}}}$

The R variables used in the preceding paragraphs are just convenient,compact expressions for the quantities shown directly above thisparagraph.

To summarize:

$\begin{matrix}{\frac{\partial F}{\partial\beta_{j}} = \left\{ \begin{matrix}{{P_{1}\beta_{1}} + {P_{2}\beta_{2}} + P_{0}} & {{{for}\mspace{14mu} j} = 1} \\{{Q_{0}(j)} + {{Q_{1}(j)}\beta_{j - 1}} + {{Q_{2}(j)}\beta_{j}} + {{Q_{3}(j)}\beta_{j + 1}}} & {{{for}\mspace{14mu} 1} < j \leq N} \\{{R_{1}\beta_{N}} + {R_{2}\beta_{N + 1}} + R_{0}} & {{{for}\mspace{14mu} j} = {N + 1}}\end{matrix} \right.} & (23)\end{matrix}$

1.3. Unity C.D.F. Constraint

From (4), we have:

$\begin{matrix}\begin{matrix}{{\int\limits_{0}^{A}{{g(x)}{x}}} = {\int\limits_{0}^{A}{\sum\limits_{i = 1}^{N}{{U_{i}(x)}{x}}}}} \\{= {\sum\limits_{i = 1}^{N}{\int\limits_{0}^{A}{{U_{i}(x)}{x}}}}} \\{= {{\sum\limits_{i = 1}^{N}{\int\limits_{a_{i}}^{a_{i + 1}}{\left( \frac{\beta_{i + 1} - \beta_{i}}{a_{i + 1} - a_{i}} \right)\left( {x - a_{i}} \right)}}} + {\beta_{i}{x}}}} \\{= {{\sum\limits_{i = 1}^{N}{\int\limits_{a_{i}}^{a_{i + 1}}{\left( \frac{\beta_{i + 1} - \beta_{i}}{a_{i + 1} - a_{i}} \right)x}}} + {\left\lbrack {\beta_{i} - \frac{a_{i}\left( {\beta_{i + 1} - \beta_{i}} \right)}{a_{i + 1} - a_{i}}} \right\rbrack {x}}}}\end{matrix} & (24)\end{matrix}$

Let:

${m_{i} = \frac{\beta_{i + 1} - \beta_{i}}{a_{i + 1} - a_{i}}},{b_{i} = {\beta_{i} - {\frac{a_{i}\left( {\beta_{i + 1} - \beta_{i}} \right)}{a_{i + 1} - a_{i\;}}.}}}$

Then (24) becomes:

$\begin{matrix}\begin{matrix}{{\int\limits_{0}^{A}{{g(x)}{x}}} = {{\sum\limits_{i = 1}^{N}{\int\limits_{a_{i}}^{a_{i + 1}}{m_{i}x}}} + {b_{i}{x}}}} \\{= {{{\sum\limits_{i = 1}^{N}{m_{i}\frac{x^{2}}{2}}} + {b_{i}x}}_{a_{i}}^{a_{i + 1}}}} \\{= {{\sum\limits_{i = 1}^{N}{m_{i}\left( {\frac{a_{i + 1}^{2}}{2} - \frac{a_{i}^{2}}{2}} \right)}} + {b_{i}\left( {a_{i + 1} - a_{i\;}} \right)}}} \\{= {{\sum\limits_{i = 1}^{N}{\left( \frac{\beta_{i + 1} - \beta_{i}}{a_{i + 1} - a_{i}} \right)\left( {\frac{a_{i + 1}^{2}}{2} - \frac{a_{i}^{2}}{2}} \right)}} +}} \\{{\left\lbrack {\beta_{i} - \frac{a_{i}\left( {\beta_{i + 1} - \beta_{i}} \right)}{a_{i + 1} - a_{i}}} \right\rbrack \left( {a_{i + 1} - a_{i}} \right)}}\end{matrix} & (25)\end{matrix}$

After further simplification, (25) becomes:

${\int\limits_{0}^{A}{{g(x)}{x}}} = {{\sum\limits_{i = 1}^{N}{T_{i}^{1}\beta_{i}}} + {T_{i}^{2}\beta_{i + 1}}}$

Where:

$\begin{matrix}{T_{i}^{1} = {a_{i + 1} - \frac{a_{i + 1}^{2} - a_{i}^{2}}{2\left( {a_{i + 1} - a_{i}} \right)}}} & (26) \\{T_{i}^{2} = {\frac{a_{i + 1}^{2} - a_{i}^{2}}{2\left( {a_{i + 1} - a_{i}} \right)} - a_{i}}} & (27)\end{matrix}$

Define:

$\begin{matrix}{{g_{1}\left( {\beta_{1},{\ldots \mspace{14mu} B_{N + 1}}} \right)} = {{\sum\limits_{i = 1}^{N}{T_{i}^{1}\beta_{i}}} + {T_{i}^{2}\beta_{i + 1}}}} & (28)\end{matrix}$

Then the unity c.d.f. constraint, (4), is expressed as

g ₁(β₁, . . . β_(N+1))=1  (29)

1.4. Constant Power Constraint

Consider (5):

$\begin{matrix}{{\int_{0}^{A}{x^{2}{g(x)}\ {x}}} = {\int_{0}^{A}{x^{2}{\sum\limits_{i = 1}^{N}{{U_{i\;}(x)}\ {x}}}}}} \\{= {\sum\limits_{i = 1}^{N}{\int_{0}^{A}{x^{2}{U_{i\;}(x)}\ {x}}}}} \\{= {\sum\limits_{i = 1}^{N}{\int_{a_{i}}^{a_{i + 1}}{x^{2}{U_{i\;}(x)}\ {x}}}}} \\{= {\sum\limits_{i = 1}^{N}{\int_{a_{i}}^{a_{i + 1}}{{x^{2}\left( {{m_{i}x} + b_{i}} \right)}{x}}}}} \\{= {{{\sum\limits_{i = 1}^{N}{m_{i}\frac{x^{4}}{4}}} + {b_{i}\frac{x^{3}}{3}}}|_{a_{i}}^{a_{i + 1}}}} \\{= {{\sum\limits_{i = 1}^{N}{m_{i}\left( {\frac{a_{i + 1}^{4}}{4} - \frac{a_{i}^{4}}{4}} \right)}} + {b_{i}\left( {\frac{a_{i + 1}^{3}}{3} - \frac{a_{i}^{3}}{3}} \right)}}}\end{matrix}$

This equation comes directly from equation 5, where we replaced the g(x)with the assumed form shown earlier. After much simplification, thisequation becomes:

${\int_{0}^{A}{x^{2}{g(x)}{x}}} = {{\sum\limits_{i = 1}^{N}{W_{i}^{1}\beta_{i}}} + {W_{i}^{2}\beta_{i + 1}}}$

Where:

$W_{i}^{1} = {\frac{a_{i + 1}^{3} - a_{i}^{3}}{3} - \frac{a_{i + 1}^{4} - a_{i}^{4}}{4\left( {a_{i + 1} - a_{i}} \right)} + \frac{a_{i}\left( {a_{i + 1}^{3} - a_{i}^{3}} \right)}{3\left( {a_{i + 1} - a_{i}} \right)}}$$W_{i}^{2} = {\frac{a_{i + 1}^{4} - a_{i}^{4}}{4\left( {a_{i + 1} - a_{i}} \right)} - \frac{a_{i}\left( {a_{i + 1}^{3} - a_{i}^{3}} \right)}{3\left( {a_{i + 1} - a_{i}} \right)}}$

Define:

$\begin{matrix}{{g_{2}\left( {B_{1},{\ldots \mspace{14mu} \beta_{N + 1}}} \right)} = {{\sum\limits_{i = 1}^{N}{W_{i}^{1}\beta_{i}}} + {W_{i}^{2}\beta_{i + 1}}}} & (30)\end{matrix}$

Then, the constant power constraint (5) becomes:

g ₂(β₁, . . . β_(N+1))=σ².

1.5. Lagrangian Definition

When we take the partial derivatives of the Lagrangian and set them tozero to find the stationary points, we get a system of equations whichincludes the original system of equations (i.e., for F), plus the twoconstraint equations. Hence, the Lagrangian contains or includes orencompasses the constraints along with the original system, allowing theLagrangian, i.e. the mathematical function that summarizes the dynamicsof the system, to be defined as:

λ(β₁, . . . β_(N+1),β₁,β₂)=F(β₁, . . . β_(N+1))+λ₁ [g ₁(β₁, . . .β_(N+1))−1]+λ₂ [g ₂(β₁, . . . β_(N+1))−σ²]

Solve the set of equations:

$\begin{matrix}{{\frac{\partial\Lambda}{\partial\beta_{j}} = 0}{\frac{\partial\Lambda}{\partial\lambda_{1}} = 0}{\frac{\partial\Lambda}{\partial\lambda_{2}} = 0.}} & (31)\end{matrix}$

We have:

$\begin{matrix}{{\frac{\partial\Lambda}{\partial\beta_{j}} = {\frac{\partial F}{\partial\beta_{j}} + {\lambda_{1}\frac{\partial g_{1}}{\partial\beta_{j}}} + {\lambda_{2}\frac{\partial g_{2}}{\partial\beta_{j}}}}}{\frac{\partial\Lambda}{\partial\lambda_{1}} = {{g_{1}\left( {\beta_{1},{\ldots \mspace{14mu} B_{N + 1}}} \right)} - 1}}{\frac{\partial\Lambda}{\partial\lambda_{2}} = {{g_{2}\left( {\beta_{1},{\ldots \mspace{14mu} B_{N + 1}}} \right)} - {\sigma^{2}.}}}} & (32)\end{matrix}$

From (32), we see that we need the partial derivatives of g₁,g₂. We setthe partial derivatives to zero, then solve for the betas to give valuesthat under certain conditions represent the maximum or minimum of thefunction. We do this for the Lagrangian, because we want to include theconstraints also, so we solve for beta parameter values which maximizeor minimize the system and which also satisfy the constraints. From(28), we have, for j=1:

$\begin{matrix}{\frac{\partial g_{1}}{\partial\beta_{j}} = \left\{ \begin{matrix}T_{1}^{1} & {{{for}\mspace{14mu} j} = 1} \\{T_{j - 1}^{2} + T_{j}^{1}} & {{{for}\mspace{14mu} 1} < j \leq N} \\T_{N}^{2} & {{{for}\mspace{14mu} j} = {N + 1}}\end{matrix} \right.} & (33)\end{matrix}$

From (30), we have:

$\begin{matrix}{\frac{\partial g_{2}}{\partial\beta_{j}} = \left\{ \begin{matrix}W_{1}^{1} & {{{{for}\mspace{14mu} j} = 1}\;} \\{W_{j - 1}^{2} + W_{j}^{1}} & {{{for}\mspace{14mu} 1} < j \leq N} \\W_{N}^{2} & {{{for}\mspace{14mu} j} = {N + 1}}\end{matrix} \right.} & (34)\end{matrix}$

Combining all the partial derivative expressions from (23), (33) and(34) gives:

$\begin{matrix}{\frac{\partial\Lambda}{\partial\beta_{j}} = \left\{ \begin{matrix}{\left( {P_{0} + {P_{1}\beta_{1}} + {P_{2}\beta_{2}}} \right) + {\lambda_{1}T_{1}^{1}} + {\lambda_{2}W_{1}^{1}}} & {{{for}\mspace{14mu} j} = 1} \\\begin{matrix}{\left\lbrack {{Q_{0}(j)} + {{Q_{1}(j)}\beta_{j - 1}} + {{Q_{2}(j)}\beta_{j}} + {{Q_{3}(j)}\beta_{j + 1}}} \right\rbrack +} \\{{\lambda_{1}\left( {T_{j - 1}^{2} + T_{j}^{1}} \right)} + {\lambda_{2}\left( {W_{j - 1}^{2} + W_{j}^{1}} \right)}}\end{matrix} & {{{for}\mspace{14mu} 1} < j \leq N} \\{\left( {R_{0\;} + {R_{1}\beta_{N}} + {R_{2}\beta_{N + 1}}} \right) + {\lambda_{1}T_{N}^{2}} + {\lambda_{2}W_{N}^{2}}} & {{{for}\mspace{14mu} j} = {N + 1}}\end{matrix} \right.} & (35)\end{matrix}$

From (35) and the system (31), we have the system of equations:

  (T₁¹)λ + (W₁¹)λ₂ + (P₁)β₁ + (P₂)β₂ = −P₀     ⋮(T_(j − 1)² + T_(j)¹)λ₁ + (W_(j − 1)² + W_(j)¹)λ₂ + (Q₁(j))β_(j − 1) + (Q₂(j))β_(j) + (Q₃(j))β_(j + 1) = −Q₀(j)     ⋮  (T_(N)²)λ₁ + (W_(N)²)λ₂ + (R₁)β_(N) + (R₂)β_(N + 1) = −R₀(0)λ₁ + (0)λ₂ + T₁¹β₁ + (T₁² + T₂¹)β₂ + … + (T_(N − 1)² + T_(N)¹)β_(N) + (T_(N)²)β_(N + 1) = 1(0)λ₁ + (0)λ₂ + W₁¹β₁ + (W₁² + W₂¹)β₂ + … + (W_(N − 1)² + W_(N)¹)β_(N) + (W_(N)²)β_(N + 1) = σ².

The vertical three dots, as used above, signify that you continue onwith the same pattern for the next value of the index j. So, the firstequation gives the expression for j=1, then you continue with the samepattern for the remaining j values until the next expression. The reasonwhy the three dots are used is two-fold: first it is compact notation,and second you don't have an exact number for the maximum value of j (itis a variable) so you can't write out all of the equations, in general.

In matrix form, the system of equations becomes:

$\begin{matrix}{{\begin{pmatrix}T_{1}^{1} & W_{1}^{1} & P_{1} & P_{2} & 0 & \; & \; & \; & \; & 0 \\{T_{1}^{2} + T_{2}^{1}} & {W_{1}^{2} + W_{2}^{1}} & {Q_{1}(2)} & {Q_{2}(2)} & {Q_{3}(2)} & 0 & \; & \; & \; & 0 \\\vdots & \; & \; & \; & \; & \; & \; & \; & \; & \; \\{T_{j - 1}^{2} + T_{j}^{1}} & {W_{j - 1}^{2} + W_{j}^{1}} & 0 & \; & 0 & {Q_{1}(j)} & {Q_{2}(j)} & {Q_{3}(j)} & 0 & 0 \\\vdots & \; & \; & \; & \; & \; & \; & \; & \; & \; \\T_{N}^{2} & W_{N}^{2} & 0 & \; & \; & \; & \; & 0 & R_{1} & R_{2} \\0 & 0 & T_{1}^{1} & {T_{1}^{2} + T_{2}^{1}} & {T_{2}^{2} + T_{3}^{1}} & \ldots & \; & \; & {T_{N - 1}^{2} + T_{N}^{1}} & T_{N}^{2} \\0 & 0 & W_{1}^{1} & {W_{1}^{2} + W_{2}^{1}} & {W_{2}^{2} + W_{3}^{1}} & \ldots & \; & \; & {W_{N - 1}^{2} + W_{N}^{1}} & W_{N}^{2}\end{pmatrix}\begin{pmatrix}\lambda_{1} \\\lambda_{2} \\\; \\\beta_{j} \\\; \\\beta_{N - 1} \\\beta_{N} \\\beta_{N + 1}\end{pmatrix}} = \begin{pmatrix}{- P_{0}} \\{- {Q_{0}(2)}} \\\vdots \\{- {Q_{0}(j)}} \\\vdots \\{- R_{0}} \\1 \\\sigma^{2}\end{pmatrix}} & (36)\end{matrix}$

The matrix form can then be solved numerically.

1.6. Hessian Matrix

A Hessian Matrix, sometimes referred to as simply a Hessian, is a squarematrix of second-order partial derivatives of a scalar-valued function,or scalar field. It describes the local curvature of a function of manyvariables. This was also previously known by the term “functionaldeterminants”.

In this case, a Hessian Matrix is useful in determination of thecritical point type, which, after solving the system (36), can beobtained from the sign test on the sequence of principal minors of thebordered Hessian [23]. To calculate the Hessian, we need the secondpartial derivatives. For the Lagrangian system described above, theHessian becomes the bordered Hessian, calculated as:

$H = \begin{pmatrix}0 & 0 & \frac{\partial g_{1}}{\partial\beta_{1}} & \ldots & \frac{\partial g_{1}}{\partial\beta_{j}} & \ldots & \frac{\partial g_{1}}{\partial\beta_{N + 1}} \\0 & 0 & \frac{\partial g_{2}}{\partial\beta_{1}} & \ldots & \frac{\partial g_{2}}{\partial\beta_{j}} & \ldots & \frac{\partial g_{2}}{\partial\beta_{N + 1}} \\\; & \; & \ddots & \; & \; & \; & \; \\\; & \; & \; & \ddots & \; & {\frac{\partial^{2}F}{{\partial\beta_{j}}{\partial\beta_{i}}} + {\lambda_{1}\frac{\partial^{2}g_{1}}{{\partial\beta_{j}}{\partial\beta_{i}}}} + {\lambda_{2}\frac{\partial^{2}g_{2}}{{\partial\beta_{j}}{\partial\beta_{i}}}}} & \; \\\; & \; & \; & \; & \ddots & \; & \;\end{pmatrix}$

From (33) and (34), we see that:

$\frac{\partial^{2}g_{1}}{{\partial\beta_{j}}{\partial\beta_{i}}} = {\frac{\partial^{2}g_{2}}{{\partial\beta_{j}}{\partial\beta_{i}}} = 0}$

From (23) we see that:

${\frac{\partial^{2}F}{\partial\beta_{1}^{2}} = P_{1}},{\frac{\partial^{2}F}{{\partial\beta_{2}}{\partial\beta_{1}}} = P_{2}},{\frac{\partial^{2}F}{{\partial\beta_{j}}{\partial\beta_{1}}} = {{0\mspace{14mu} {for}\mspace{14mu} j} > 2}}$${\frac{\partial^{2}F}{\partial\beta_{j}^{2}} = {Q_{2}(j)}},{\frac{\partial^{2}F}{{\partial\beta_{j - 1}}{\partial\beta_{j}}} = {Q_{1}(j)}},{\frac{\partial^{2}F}{{\partial\beta_{j + 1}}{\partial\beta_{j}}} = {Q_{3}(j)}},{\frac{\partial^{2}F}{{\partial\beta_{k}}{\partial\beta_{j}}} = {{0\mspace{14mu} {for}\mspace{14mu} k} \neq \left\{ {{j - 1},j,{j + 1}} \right\}}},{\frac{\partial^{2}F}{\partial\beta_{N + 1}^{2}} = R_{2}},{\frac{\partial^{2}F}{{\partial\beta_{N}}{\partial\beta_{N + 1}}} = R_{1}},{\frac{\partial^{2}F}{{\partial\beta_{j}}{\partial\beta_{N + 1}}} = {{0\mspace{14mu} {for}\mspace{14mu} j} < N}}$

Performing the sign test numerically shows that the solution to (36) isa minimum.

1.7 Numerical Solutions and Analysis

We now can solve system (36) numerically. Throughout, we take σ=√{squareroot over (2)}/2 to generate a normalized complex signal with unitpower. FIGS. 2 through 9 contain plots of the optimal solution comparedagainst the baseline Rayleigh distribution, as a function of cutoffamplitude value A. Values of N range from 8 to 11, depending on A. Abovethese N values, numerical problems occur causing the determinant of thesystem to be close to zero and thus precluding numerical solution.

Two observations are clear from the plots in FIGS. 2 through 9. Thefirst observation is that as the value of A increases the solutionapproaches the baseline Rayleigh distribution. This is most apparentwhen comparing the solution prior to the tail of the Rayleighdistribution, since most of the solutions appear to match the Rayleighdistribution reasonably well before the tail. On the tail, the solutiontends to curve upward to satisfy the unity c.d.f. and constant powerconstraints. This behavior suggests that closely matching the Rayleighdistribution when possible, is a desirable solution property.

The second observation is that for almost all of the solutions, thefirst ordinate, i.e. the y-coordinate, representing the distance from apoint to the horizontal or x-axis measured parallel to the vertical or yaxis, in this case β₁, is negative and thus a portion of the solutionhas negative values. This condition violates the non-negativityconstraint on probability density functions.

There are two methods for overcoming the partial negativity of theoptimal solution. The first method increases the first abscissa, theperpendicular distance of a point from the vertical axis, a₁=0 value toa small but positive value, specifically until β₁ in the solutionbecomes positive. In this case, the solution consists of the piecewisesolution for x≧a₁, and the up-front portion of the Rayleigh distributionfor 0≦x<a₁. The perturbation of a₁ has generally been observed to besmall.

As an example, FIG. 10 contains the solution for A=1.4. In thatinstance, a value of a₁=0.0082 is sufficient to make β₁>0. The solutionshown in FIG. 10 is negligibly different from the optimal solution shownin FIG. 5.

The second method for overcoming the negativity of the solution is touse the Rayleigh distribution as the front segment of the distribution,and then use a perturbed version of the back-end of the optimalsolution. In this approach, the tail of the optimal solution must beperturbed to meet the new unity c.d.f. and constant power constraints,which are generated by replacing the front-end of the optimal solutionwith the Rayleigh distribution. This perturbation approach is discussedin the next section.

1.8. Perturbations from the Optimal Solution

Solutions of the system (36) tend to give a negative value for the firstordinate value β₁. Therefore, to determine a solution for which thenon-negativity condition (6) holds, we perturb the optimal solution. Wechoose a_(k) starting at the distribution tail to perturb the solution.For xε[0,a₁] we use the original Rayleigh distribution ƒ_(R)(x), and weperturb the ordinate values for a_(j) with j>k for some k<N as depictedin FIG. 11.

We now define a perturbed function as:

$\begin{matrix}{{g_{ɛ}(x)} = {\sum\limits_{i = k}^{N}{U_{i}\left( {x,ɛ_{i}} \right)}}} & (37)\end{matrix}$

We define the perturbed function this way, because we intend to changethe beta values by small amounts, i.e. by the epsilon values.

For some 1<k<N, where the U_(i)(x,ε_(i)) are perturbed,compactly-supported linear segments:

$\begin{matrix}{{U_{i}\left( {x,ɛ_{i}} \right)} = \left\{ \begin{matrix}{{\left( \frac{\left( {\beta_{i + 1} + ɛ_{i + 1}} \right) - \left( {\beta_{i} + ɛ_{i}} \right)}{a_{i + 1} - a_{i}} \right)\left( {x - a_{i}} \right)} + \left( {\beta_{i} + ɛ_{i}} \right)} & {{{for}\mspace{14mu} x} \in \left\lbrack {a_{i},a_{i + 1}} \right\rbrack} \\0 & {otherwise}\end{matrix} \right.} & (38)\end{matrix}$

Here, k represents the index value of the a value immediately before thea values corresponding to the beta values that we are going to change,see FIG. 11. We pick an index k, past which we adjust the betaparameters.

We then take ε_(k)=0 so that (a_(k), β_(k)+ε_(k))=(a_(k),β_(k)) thusleaving the point at index k fixed, maintaining a continuous function.Unity c.d.f and constant power constraints must still be satisfied forthe perturbed solution; these are examined next.

1.8.1 Unity C.D.F. Constraint

From the unity c.d.f. condition on g_(ε)(x), we have:

$\begin{matrix}{1 = {\int_{0}^{A}{{g_{ɛ}(x)}\ {x}}}} \\{= {{\int_{0}^{a_{k}}{{f_{R}(x)}\ {x}}} + {\int_{a_{k}}^{A}{{g_{ɛ}(x)}\ {x}}}}}\end{matrix}$

which gives:

$\begin{matrix}{{\int_{a_{k}}^{A}{{g_{ɛ}(x)}\ {x}}} = {1 - {\int_{0}^{a_{k}}{{f_{R}(x)}\ {x}}}}} & (39)\end{matrix}$

From integration of the Rayleigh distribution in (16), we obtain:

$\begin{matrix}{{\int_{a_{k}}^{A}{{g_{ɛ}(x)}\ {x}}} = {1 - \left\lbrack {1 - ^{- \frac{a_{k}^{2}}{\sigma^{2}}}} \right\rbrack}} \\{= ^{- \frac{a_{k}^{2}}{\sigma^{2}}}}\end{matrix}$

Now from (37):

$\begin{matrix}{{\int_{a_{k}}^{A}{{g_{ɛ}(x)}\ {x}}} = {\int_{a_{k}}^{A}{\sum\limits_{i = k}^{N}{{U_{i}\left( {x,ɛ_{i}} \right)}\ {x}}}}} \\{= {\sum\limits_{i = k}^{N}{\int_{a_{k}}^{A}{{U_{i}\left( {x,ɛ_{i}} \right)}\ {x}}}}} \\{= {\sum\limits_{i = k}^{N}{\int_{a_{i}}^{a_{i + 1}}{{U_{i}\left( {x,ɛ_{i}} \right)}\ {x}}}}} \\{= {{\sum\limits_{i = k}^{N}{\int_{a_{i}}^{a_{i + 1}}{m_{i}x}}} + {b_{i}\ {x}}}}\end{matrix}$

Where:

$m_{i} = \frac{\left( {\beta_{i + 1} - \beta_{i}} \right) + \left( {ɛ_{i + 1} - ɛ_{i}} \right)}{a_{i + 1} - a_{i}}$$b_{i} = {\beta_{i} + ɛ_{i} - {\frac{a_{i}\left\lbrack {\left( {\beta_{i + 1} - \beta_{i}} \right) + \left( {ɛ_{i + 1} - ɛ_{i}} \right)} \right\rbrack}{a_{i + 1} - a_{i}}.}}$

Thus:

$\begin{matrix}{{{{\sum\limits_{i = k}^{N}{\int_{a_{i}}^{a_{i + 1}}{m_{i}x}}} + {b_{i}\ {x}}} = {{\sum\limits_{i = k}^{N}{m_{i}\frac{x^{2}}{2}}} + {b_{i}x}}}}_{a_{i}}^{a_{i + 1}} \\{= {{\sum\limits_{i = k}^{N}{m_{i}\left( {\frac{a_{i + 1}^{2}}{2} - \frac{a_{i}^{2}}{2}} \right)}} + {b_{i}\left( {a_{i + 1} - a_{i}} \right)}}}\end{matrix}$

Which, after much simplification, gives:

$\begin{matrix}{{\int_{a_{k}}^{A}{{g_{ɛ}(x)}\ {x}}} = {{\sum\limits_{i = k}^{N}{T_{i}^{1}\beta_{i}}} + {T_{i}^{2}\beta_{i + 1}} + {ɛ_{i}\left( \frac{a_{i + 1} - a_{i}}{2} \right)} + {ɛ_{i + 1}\left( \frac{a_{i + 1} - a_{i}}{2} \right)}}} & (40)\end{matrix}$

for T_(i) ¹ and T_(i) ² defined in (26) and (27).

Hence from (39) and (40) we get:

$\begin{matrix}{{{\sum\limits_{i = k}^{N}{ɛ_{i}\left( \frac{a_{i + 1} - a_{i}}{2} \right)}} + {ɛ_{i + 1}\left( \frac{a_{i + 1} - a_{i}}{2} \right)}} = {1 - {\int_{0}^{a_{k}}{{f_{R}(x)}\ {x}}} - {\sum\limits_{i = k}^{N}{T_{i}^{1}\beta_{i}}} + {T_{i}^{2}\beta_{i + 1}}}} & (41)\end{matrix}$

1.8.2 Constant Power Constraint

In general, we want to keep the average power level the same as for theoriginal, uncompanded OFDM signal, while simultaneously reducing thepeak power, so as to lower the PAPR value of the companded OFDM signalrelative to the nominal OFDM signal, so a constant power constraint isused.

From the constant power constraint on g_(ε)(x), we have:

$\begin{matrix}{\sigma^{2} = {\int_{0}^{A}{x^{2}{g_{ɛ}(x)}\ {x}}}} \\{= {{\int_{0}^{a_{k}}{x^{2}{f_{R}(x)}\ {x}}} + {\int_{a_{k}}^{A}{x^{2}{g_{ɛ}(x)}\ {x}}}}}\end{matrix}$

which gives:

$\begin{matrix}{{\int_{a_{k}}^{A}{x^{2}{g_{ɛ}(x)}\ {x}}} = {\sigma^{2} - {\int_{0}^{a_{k}}{x^{2}{f_{R}(x)}\ {x}}}}} & (42)\end{matrix}$

Using integration by parts, we get:

${\int_{c}^{d}{x^{2}{f_{R}(x)}\ {x}}} = {{\left( {c^{2} + \sigma^{2}} \right)^{- \frac{c^{2}}{\sigma^{2}}}} - {\left( {d^{2} + \sigma^{2}} \right)^{\frac{d^{2}}{\sigma^{2}}}}}$

and so:

${\int_{0}^{a_{k}}{x^{2}{f_{R}(x)}\ {x}}} = {\sigma^{2} - {\left( {a_{k}^{2} + \sigma^{2}} \right)^{- \frac{a_{k}^{2}}{\sigma^{2}}}}}$

from which follows, using (42):

${\int_{a_{k}}^{A}{x^{2}{g_{ɛ}(x)}\ {x}}} = {\left( {a_{k}^{2} + \sigma^{2}} \right)^{- \frac{a_{k}^{2}}{\sigma^{2}}}}$

Now:

$\begin{matrix}{{\int_{a_{k}}^{A}{x^{2}{g_{ɛ}(x)}\ {x}}} = {\int_{a_{k}}^{A}{x^{2}{\sum\limits_{i = k}^{N}\ {{U_{i}\left( {x,ɛ_{i}} \right)}{x}}}}}} \\{= {\sum\limits_{i = k}^{N}{\int_{a_{k}}^{A}{x^{2}\ {U_{i}\left( {x,ɛ_{i}} \right)}{x}}}}} \\{= {\sum\limits_{i = k}^{N}{\int_{a_{i}}^{a_{i + 1}}{x^{2}\ {U_{i}\left( {x,ɛ_{i}} \right)}{x}}}}} \\{= {\sum\limits_{i = k}^{N}{\int_{a_{i}}^{a_{i + 1}}{{x^{2}\left( {{m_{i}x} + b_{i}} \right)}\ {x}}}}}\end{matrix}$

Which, after much manipulation, becomes:

$\begin{matrix}{{\int_{a_{k}}^{A}{x^{2}{g_{ɛ}(x)}\ {x}}} = {{\sum\limits_{i = k}^{N}{W_{i}^{1}\beta_{i}}} + {W_{i}^{2}\beta_{i + 1}} + {ɛ_{i}\eta_{i}^{1}} + {ɛ_{i + 1}\eta_{i}^{2}}}} & (43)\end{matrix}$

Where:

$\eta_{i}^{1} = {{\frac{a_{i + 1}^{3} - a_{i}^{3}}{3} - \frac{a_{i + 1}^{4} - a_{i}^{4}}{4\left( {a_{i + 1} - a_{i}} \right)} + \frac{a_{i}\left( {a_{i + 1}^{3} - a_{i}^{3}} \right)}{3\left( {a_{i + 1} - a_{i}} \right)}} = W_{i}^{1}}$$\eta_{i}^{2} = {{\frac{a_{i + 1}^{4} - a_{i}^{4}}{4\left( {a_{i + 1} - a_{i}} \right)} - \frac{a_{i}\left( {a_{i + 1}^{3} - a_{i}^{3}} \right)}{3\left( {a_{i + 1} - a_{i}} \right)}} = {W_{i}^{2}.}}$

From (43) using (42), we get:

$\begin{matrix}{{{\sum\limits_{i = k}^{N}{ɛ_{i}\eta_{i}^{1}}} + {ɛ_{i + 1}\eta_{i}^{2}}} = {\sigma^{2} - {\int_{0}^{a_{k}}{x^{2}{f_{R}(x)}\ {x}}} - {\sum\limits_{i = k}^{N}{W_{i}^{1}\beta_{i}}} + {W_{i}^{2}\beta_{i + 1}}}} & (44)\end{matrix}$

Combining (41) and (44), we get the under-constrained system ofequations:

$\begin{matrix}{{{\sum\limits_{i = k}^{N}{\xi_{i}ɛ_{i}}} + {\xi_{i}ɛ_{i + 1}}} = C_{1}} & (45) \\{{{\sum\limits_{i = k}^{N}{\eta_{i}^{1}ɛ_{i}}} + {\eta_{i}^{2}ɛ_{i + 1}}} = C_{2}} & (46)\end{matrix}$

Where:

$\xi_{i} = \frac{a_{i + 1} - a_{i}}{2}$$C_{1} = {1 - {\int_{0}^{a_{k}}{f_{R}\; (x)\ {x}}} - {\sum\limits_{i = k}^{N}{T_{i}^{1}\beta_{i}}} + {T_{i}^{2}\beta_{i + 1}}}$$C_{2} = {\sigma^{2} - {\int_{0}^{a_{k}}{x^{2}f_{R}\; (x)\ {x}}} - {\sum\limits_{i = k}^{N}{W_{i}^{1}\beta_{i}}} + {W_{i}^{2}{\beta_{i + 1}.}}}$

Next, the optimal values of ε_(i) are found.

1.8.3 Optimal Perturbation Parameters

Our goal here is to find the minimal perturbation g_(ε)(x) away fromg(x). To do so, we'll choose parameters (ε_(k+1), . . . ε_(N)) tominimize:

$\begin{matrix}{{F_{ɛ}\left( {{ɛ_{{k + 1},}\mspace{11mu} \ldots}\mspace{14mu},ɛ_{N + 1}} \right)} = {\sum\limits_{i = {k + 1}}^{N + 1}ɛ_{i}^{2}}} & (47)\end{matrix}$

First we rewrite (45) and (46) as:

$\begin{matrix}\begin{matrix}{C_{1} = {{\sum\limits_{i = k}^{N}{\xi_{i}ɛ_{i}}} + {\xi_{i}ɛ_{i + 1}}}} \\{= {{\xi_{k}ɛ_{k}} + \left( {{\sum\limits_{i = {k + 1}}^{N - 1}{\xi_{i - 1}ɛ_{i}}} + {\xi_{i}ɛ_{i}}} \right) + \left( {{\xi_{N - 1}ɛ_{N}} + {\xi_{N}ɛ_{N}}} \right) + {\xi_{N}ɛ_{N + 1}}}} \\{= {{\xi_{k}ɛ_{k}} + \left( {\sum\limits_{i = {k + 1}}^{N - 1}{\left( {\xi_{i - 1} + \xi_{i}} \right)ɛ_{i}}} \right) + {\left( {\xi_{N - 1} + \xi_{N}} \right)ɛ_{N}} + {\xi_{N}ɛ_{N + 1}}}} \\{= {\left( {\sum\limits_{i = {k + 1}}^{N - 1}{\left( {\xi_{i - 1} + \xi_{i}} \right)ɛ_{i}}} \right) + {\left( {\xi_{N - 1} + \xi_{N}} \right)ɛ_{N}} + {\xi_{N}ɛ_{N + 1}}}}\end{matrix} & (48)\end{matrix}$

because, by definition, we took ε_(k)=0. Similarly from (46):

$\begin{matrix}\begin{matrix}{C_{2} = {{\sum\limits_{i = k}^{N}{\eta_{i}^{1}ɛ_{i}}} + {\eta_{i}^{2}ɛ_{i + 1}}}} \\{= {{\eta_{k}^{1}ɛ_{k}} + \left( {{\sum\limits_{i = {k + 1}}^{N - 1}{\eta_{i - 1}^{2}ɛ_{i}}} + {\eta_{i}^{1}ɛ_{i}}} \right) + \left( {{\eta_{N - 1}^{2}ɛ_{N}} + {\eta_{N}^{1}ɛ_{N}}} \right) + {\eta_{N}^{2}ɛ_{N + 1}}}} \\{= {{\eta_{k}^{1}ɛ_{k}} + \left( {\sum\limits_{i = {k + 1}}^{N - 1}{\left( {\eta_{i - 1}^{2} + \eta_{i}^{1}} \right)ɛ_{i}}} \right) + {\left( {\eta_{N - 1}^{2} + \eta_{N}^{1}} \right)ɛ_{N}} + {\eta_{N}^{2}ɛ_{N + 1}}}} \\{= {\left( {\sum\limits_{i = {k + 1}}^{N - 1}{\left( {\eta_{i - 1}^{2} + \eta_{i}^{1}} \right)ɛ_{i}}} \right) + {\left( {\eta_{N - 1}^{2} + \eta_{N}^{1}} \right)ɛ_{N}} + {\eta_{N}^{2}ɛ_{N + 1}}}}\end{matrix} & (49)\end{matrix}$

Now we define functions:

$\begin{matrix}{{h_{1}\left( {ɛ_{k + 1},\ldots \mspace{14mu},ɛ_{N + 1}} \right)} = {\left( {\sum\limits_{i = {k + 1}}^{N - 1}{\left( {\xi_{i - 1} + \xi_{i}} \right)ɛ_{i}}} \right) + {\left( {\xi_{N - 1} + \xi_{N}} \right)ɛ_{N}} + {\xi_{N}ɛ_{N + 1}}}} & (50) \\{{h_{2}\left( {ɛ_{k + 1},\ldots \mspace{14mu},ɛ_{N + 1}} \right)} = {\left( {\sum\limits_{i = {k + 1}}^{N - 1}{\left( {\eta_{i - 1}^{2} + \eta_{i}^{1}} \right)ɛ_{i}}} \right) + {\left( {\eta_{N - 1}^{2} + \eta_{N}^{1}} \right)ɛ_{N}} + {\eta_{N}^{2}ɛ_{N + 1}}}} & (51)\end{matrix}$

Then from (48) and (49) we get two constraint equations:

h ₁(ε_(k+1), . . . ε_(N+1))=C ₁

h ₂(ε_(k+1), . . . ε_(N+1))=C ₂

Next, we define the Lagrangian:

Λ_(h)(ε_(k+1), . . . ,ε_(N+1),λ₁,λ₂)=F _(ε)(ε_(k+1), . . . ,ε_(N+1))+λ₁[h ₁(ε_(k+1), . . . ,ε_(N+1))C ₁]+λ₂ [h ₂(ε_(k+1), . . . ε_(N+1))−C ₂]

Now, we solve the set of equations:

$\begin{matrix}{{\frac{\partial\Lambda_{h}}{\partial ɛ_{j}} = 0}{\frac{\partial\Lambda_{h}}{\partial\lambda_{1}} = 0}{\frac{\partial\Lambda_{h}}{\partial\lambda_{2}} = 0.}} & (52)\end{matrix}$

We have:

$\frac{\partial\Lambda_{h}}{\partial ɛ_{j}} = {\frac{\partial F_{ɛ}}{\partial ɛ_{j}} + {\lambda_{1}\frac{\partial h_{1}}{\partial ɛ_{j}}} + {\lambda_{2}\frac{\partial h_{2}}{\partial ɛ_{j}}}}$$\frac{\partial\Lambda_{h}}{\partial\lambda_{1}} = {{h_{1}\left( {ɛ_{k + 1},\ldots \mspace{14mu},ɛ_{N + 1}} \right)} - C_{1}}$$\frac{\partial\Lambda_{h}}{\partial\lambda_{2}} = {{h_{2}\left( {ɛ_{k + 1},\ldots \mspace{14mu},ɛ_{N + 1}} \right)} - C_{2}}$

From (47):

$\frac{\partial F_{ɛ}}{\partial ɛ_{j}} = {2ɛ_{j}}$

From (50) and (51):

$\frac{\partial h_{1}}{\partial ɛ_{i}} = \left\{ {\begin{matrix}{\xi_{i - 1} + \xi_{i}} & {{{{if}\mspace{14mu} k} + 1} \leq i \leq N} \\\xi_{N} & {{{if}\mspace{14mu} i} = {N + 1}}\end{matrix},{\frac{\partial h_{2}}{\partial ɛ_{i}} = \left\{ \begin{matrix}{\eta_{i - 1}^{2} + \eta_{i}^{1}} & {{{{if}\mspace{14mu} k} + 1} \leq i \leq N} \\\eta_{N}^{2} & {{{if}\mspace{14mu} i} = {N + 1}}\end{matrix} \right.}} \right.$

Thus:

$\frac{\partial\Lambda_{j}}{\partial ɛ_{j}} = \left\{ \begin{matrix}{{2ɛ_{i}} + {\lambda_{1}\left( {\xi_{i - 1} + \xi_{i\;}} \right)} + {\lambda_{2}\left( {\eta_{i - 1}^{2} + \eta_{i}^{1}} \right)}} & {{{{if}\mspace{14mu} k} + 1} \leq i \leq N} \\{{2ɛ_{N + 1}} + {\lambda_{1}\left( \xi_{N} \right)} + {\lambda_{2}\left( \eta_{N}^{2} \right)}} & {{{if}\mspace{14mu} i} = {N + 1}}\end{matrix} \right.$

Hence, we need to solve the following system of equations for ε_(i), λ₁,λ₂:

  (ξ_(k) + ξ_(k + 1))λ₁ + (η_(k)² + η_(k + 1)¹)λ₂ + (2)ɛ_(k + 1) = 0     ⋮  (ξ_(i − 1) + ξ_(i))λ₁ + (η_(i − 1)² + η_(i)¹)λ₂ + (2)ɛ_(i) = 0     ⋮  (ξ_(N))λ₁ + (η_(N)²)λ₂ + (2)ɛ_(N + 1) = 0${{(0)\lambda_{1}} + {(0)\lambda_{2}} + \left( {\sum\limits_{i = {k + 1}}^{N - 1}{\left( {\xi_{i - 1} + \xi_{i}} \right)ɛ_{i\;}}} \right) + {\left( {\xi_{N - 1} + \xi_{N}} \right)ɛ_{N}} + {\xi_{N}ɛ_{N + 1}}} = {{{{C_{1}(0)}\lambda_{1}} + {(0)\lambda_{2}} + \left( {\sum\limits_{i = {k + 1}}^{N - 1}{\left( {\eta_{i - 1}^{2} + \eta_{i}^{1}} \right)ɛ_{i\;}}} \right) + {\left( {\eta_{N - 1}^{2} + \eta_{N}^{1}} \right)ɛ_{N}} + {\eta_{N}^{2}ɛ_{N + 1}}} = C_{2}}$

Or:

$\begin{matrix}{{\begin{pmatrix}\left( {\xi_{k} + \xi_{k + 1}} \right) & \left( {\eta_{k}^{2} + \eta_{k + 1}^{1}} \right) & 2 & 0 & \; & \; & \; & \; & 0 \\\left( {\xi_{k + 1} + \xi_{k + 2}} \right) & \left( {\eta_{k + 1}^{2} + \eta_{k + 2}^{1}} \right) & 0 & 2 & 0 & \; & \; & \; & 0 \\\vdots & \vdots & \; & \; & \ddots & \; & \; & \; & \; \\\; & \; & \; & \; & \; & \ddots & \; & \; & \; \\\left( {\xi_{i + 1} + \xi_{i}} \right) & \left( {\eta_{i + 1}^{2} + \eta_{i}^{1}} \right) & 0 & \; & \; & \; & 2 & 0 & 0 \\\vdots & \vdots & \; & \; & \; & \; & \; & \ddots & \; \\\xi_{N} & \eta_{N}^{2} & 0 & \; & \; & \; & \; & 0 & 2 \\0 & 0 & \left( {\xi_{k} + \xi_{k + 1}} \right) & \; & \; & \ldots & \; & \; & \xi_{N} \\0 & 0 & \left( {\eta_{k}^{2} + \eta_{k + 1}^{1}} \right) & \; & \; & \ldots & \; & \; & \eta_{N}^{2}\end{pmatrix}\begin{pmatrix}\lambda_{1} \\\lambda_{2} \\\; \\\; \\{\; ɛ_{i}} \\\; \\ɛ_{N - 1} \\ɛ_{N} \\ɛ_{N + 1}\end{pmatrix}} = \begin{pmatrix}0 \\0 \\\vdots \\\; \\\; \\\vdots \\\; \\C_{1} \\C_{2}\end{pmatrix}} & (53)\end{matrix}$

This system of equations (53) may then be solved numerically.

1.9. Perturbed Solution Numerical Results

Perturbation solutions for the system in (53) are shown for A=1.4 andkε{5,6,7,8} in FIGS. 12 through 15. In these Figures, the perturbationsolution with the Rayleigh front-end is plotted against the optimalsolution.

1.10. Compander Derivation

Next, we denote the companding function by C, which operates on theinput sequence x(n) by modifying the amplitude to produce the compandedoutput sequence y(n)=C{x(n)}. The compander is monotonic, i.e. itmaintains order of the function (i.e., either increasing or decreasing),over its domain xε[0,A], so the compander may be derived as:

C{x(n)}=sgn{x(n)}F _(|y|) ⁻¹ F _(|x|) {x(n)}  (54)

Where F_(|x|) denotes the c.d.f of the input and F_(|y|) denotes thec.d.f of the output companded signal. To derive the compander,expressions for each c.d.f in (54) are needed. The p.d.f. of thecompanded signal is given by:

${f_{y}(x)} = \left\{ \begin{matrix}{\frac{2x}{\sigma^{2}}^{- \frac{x^{2}}{\sigma^{2}}}} & {{{for}\mspace{14mu} 0} \leq x < a_{1}} \\{f_{a_{1}}(x)} & {{{for}\mspace{14mu} a_{1}} \leq x < a_{2}} \\\vdots & \vdots \\{f_{a_{i}}(x)} & {{{for}\mspace{14mu} a_{i}} \leq x < a_{i + 1}} \\\vdots & \vdots \\{f_{a_{N}}(x)} & {{{for}\mspace{14mu} a_{N}} \leq x < a_{N + 1}} \\0 & {{{for}\mspace{14mu} a_{N + 1}} \leq x}\end{matrix} \right.$

For a_(i)≦x<a_(i+1), the c.d.f. is given by:

$\begin{matrix}\begin{matrix}{{F_{y}(x)} = {\int_{0}^{x}{{f_{y}(\xi)}\ {\xi}}}} \\{= {\left\lbrack {{\int_{0}^{a_{1}}{{f_{R}(\xi)}\ {\xi}}} + {\sum\limits_{j = 1}^{i - 1}{\int_{a_{j}}^{a_{j + 1}}{{f_{a_{j}}\ (\xi)}{\xi}}}}} \right\rbrack + {\int_{a_{i}}^{x}{{f_{a_{i}}(\xi)}\ {\xi}}}}} \\{= {Z_{i - 1} + {\int_{a_{i}}^{x}{{f_{a_{i}}(\xi)}\ {\xi}}}}}\end{matrix} & \begin{matrix}(55) \\(56)\end{matrix}\end{matrix}$

Where Z_(i−1) denotes the term in brackets. Thus:

$\begin{matrix}\begin{matrix}{{{F_{y}(x)} - Z_{i - 1}} = {\int_{0}^{x}{{f_{a_{i}}(\xi)}\ {\xi}}}} \\{= {{\int_{a_{i}}^{x}{m_{i}\xi}} + {b_{i}\ {\xi}}}} \\{= {{m_{i}\frac{x^{2}}{2}} + {b_{i}x} - {m_{i}\frac{a_{i}^{2}}{2}} - {b_{i}a_{i}}}} \\{= {{\left( \frac{m_{i}}{2} \right)x^{2}} + {\left( b_{i} \right)x} + \left\lbrack {- \left( {{m_{i}\frac{a_{i}^{2}}{2}} + {b_{i}a_{i}}} \right)} \right\rbrack}} \\{= {{(\gamma)x^{2}} + {(\delta)x} + e}}\end{matrix} & (57)\end{matrix}$

with the obvious definitions for γ, δ, e, i.e. they correspond to thecoefficients of the x times x term, the x term, and the constant term.Now, γ might be negative, so divide both sides of (57) by γ to get:

$\begin{matrix}{\frac{{F_{y}(x)} - Z_{i - 1}}{\gamma} = {x^{2} + {\left( \frac{\delta}{\gamma} \right)x} + \left( \frac{e}{\gamma} \right)}} & (58)\end{matrix}$

Complete the square in (58):

$\begin{matrix}{\left( {x + \frac{\delta}{2\gamma}} \right)^{2} = {x^{2} + {\frac{2\delta}{2\gamma}x} + \frac{\delta^{2}}{4\gamma^{2}}}} \\{= {x^{2} + {\left( \frac{\delta}{\gamma} \right)x} + \left( \frac{\delta^{2}}{4\gamma^{2}} \right)}}\end{matrix}$

Which implies:

$\begin{matrix}\begin{matrix}{{\left( {x + \frac{\delta}{2\gamma}} \right)^{2} + \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)} = {x^{2} + {\left( \frac{\delta}{\gamma} \right)x} + \frac{\delta^{2}}{4\gamma^{2}} + \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)}} \\{= {x^{2} + {\left( \frac{\delta}{\gamma} \right)x} + \left( \frac{e}{\gamma} \right)}}\end{matrix} & (59)\end{matrix}$

Using (59), (58) becomes:

$\begin{matrix}{\frac{{F_{y}(x)} - Z_{i - 1}}{\gamma} = {\left. {\left( {x + \frac{\delta}{2\gamma}} \right)^{2} + \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)}\Rightarrow \left( {x + \frac{\delta}{2\gamma}} \right)^{2} \right. = {\left. {\frac{{F_{y}(x)} - Z_{i - 1}}{\gamma} - \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)}\Rightarrow {x + \frac{\delta}{\gamma^{2}}} \right. = \sqrt[ \pm ]{\frac{{F_{y}(x)} - Z_{i - 1}}{\gamma} - \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)}}}} & (60)\end{matrix}$

(The above arrows are mathematical shorthand for the phrase: “whichimplies”)

And thus:

$\begin{matrix}{x = {\sqrt[ \pm ]{\frac{{F_{y}(x)} - Z_{i - 1}}{\gamma} - \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)} - \frac{\delta}{2\gamma}}} & (61)\end{matrix}$

In (61) we need to decide whether to take the positive or negativesquare-root. The decision is made from (60) by examining the sign of

$x + {\frac{\delta}{2\gamma}.}$

If positive, the positive square-root is chosen, if negative, thenegative square-root should be chosen. Of course, if x≧A, then weconsider the sign on

$A + {\frac{\delta}{2\gamma}.}$

Finally, the compander is given by:

$\begin{matrix}{{F_{y}^{- 1}{F_{R}(x)}} = {\sqrt[ \pm ]{\frac{1 - ^{- \frac{x^{2}}{\sigma^{2}}} - Z_{i - 1}}{\gamma} - \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)} - \frac{\delta}{2\gamma}}} & (62)\end{matrix}$

for a_(i)≦x<a_(a+1). To complete the compander derivation, we need tocalculate Z_(i) From (16) and (55):

$\begin{matrix}{Z_{0} = {\int_{0}^{a_{1}}{{f_{R}(\xi)}\ {\xi}}}} \\{= {1 - ^{- \frac{a_{1}^{2}}{\sigma^{2}}}}}\end{matrix}$

For i>0, Z_(i) contains integral terms of the form:

∫_(a_(j))^(a_(j + 1))f_(a_(j))(ξ) ξ

which represents a trapezoidal area and can be shown to be:

${\int_{a_{j}}^{a_{j + 1}}{{f_{a_{j}}(\xi)}\ {\xi}}} = \frac{{\beta_{i + 1}\left( {a_{i + 1} - a_{i}} \right)} + {\beta_{i}\left( {a_{i + 1} - a_{i}} \right)}}{2}$

Hence:

$\begin{matrix}{Z_{i} = {Z_{i - 1} + \frac{{\beta_{i + 1}\left( {a_{i + 1} - a_{i}} \right)} + {\beta_{i}\left( {a_{i + 1} - a_{i}} \right)}}{2}}} & (63)\end{matrix}$

1.11 Decompander Derivation

To find the decompander, we need to invert (62):

$\begin{matrix}{y = {\left. {\sqrt[ \pm ]{\frac{1 - ^{- \frac{x^{2}}{\sigma^{2}}} - Z_{i - 1}}{\gamma} - \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)} - \frac{\delta}{2\gamma}}\Rightarrow \left( {y + \frac{\delta}{2\gamma}} \right)^{2} \right. = {\left. {\frac{1 - ^{- \frac{x^{2}}{\sigma^{2}}} - Z_{i - 1}}{\gamma} - \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)}\Rightarrow {\gamma \left( {y + \frac{\delta}{2\gamma}} \right)}^{2} \right. = \left. {1 - ^{- \frac{x^{2}}{\sigma^{2}}} - Z_{i - 1} - e + \frac{\delta^{2}}{4\gamma}}\Rightarrow {^{- \frac{x^{2}}{\sigma^{2}}} - {\gamma \left( {y + \frac{\delta}{2\gamma}} \right)}^{2} - e + \frac{\delta^{2}}{4\gamma} + 1 - Z_{i - 1}}\Rightarrow \right.}}} & (63) \\{x = {{{sgn}(x)}\sqrt{{- \sigma^{2}}{\ln \left( {{- {\gamma \left( {y + \frac{\delta}{2\gamma}} \right)}^{2}} - e + \frac{\delta^{2}}{4\gamma} + 1 - Z_{i - 1}} \right)}}}} & (64)\end{matrix}$

Note that, in (64), the sgn function refers to retaining either thepositive or negative square-root, depending on the sign of X.

Finally, we need to find the range of the compander for use in thedecompander. When a_(i)≦x<a_(i+1), then y_(MIN)≦y<y_(MAX) where:

y _(MIN)=min{y(a _(i)),y(a _(i+1))}

y _(MAX)=max{y(a _(i)),y(a _(i+1))}

and from (63):

${y\left( a_{i} \right)} = {\sqrt[ \pm ]{\frac{1 - ^{- \frac{a_{i}^{2}}{\sigma^{2}}} - Z_{i - 1}}{\gamma} - \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)} - \frac{\delta}{2\gamma}}$${y\left( a_{i + 1} \right)} = {\sqrt[ \pm ]{\frac{1 - ^{- \frac{a_{i + 1}^{2}}{\sigma^{2}}} - Z_{i - 1}}{\gamma} - \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)} - \frac{\delta}{2\gamma}}$

2. SIMULATION

Simulation results were generated for a Quadrature Phase Shift Keying(QPSK) (a type of phase shift keying) modulated OFDM signal with sixtyfour subcarriers over an additive white Gaussian noise channel. Tenthousand realizations of signal-plus-noise were generated for each noisepower.

For the Rayleigh p.d.f., σ=√{square root over (2)}/2 was chosen. Anoversampling factor of four was used so that the discrete PAPRreasonably approximates the PAPR from a continuous system as isdescribed Y. Wang, L.-H. Wang, J.-H. Ge, and B. Ai., “An EfficientNonlinear Companding Transform for Reducing PAPR of OFDM Signals”, IEEETrans. Broadcast., vol. 58, no. 4, pp.: 677-684, December 2012. PAPRreduction was measured using the complimentary cumulative distributionfunction (CCDF), which gives the probability that the PAPR is above agiven threshold value. Symbol error rates, to measure demodulationperformance, and power spectrums, to measure out-of-band powerreduction, were also generated.

Lagrange compander results were then generated and compared against twostate-of-the-art companders; the linear compander (LIN) in Y. Wang,L.-H. Wang, J.-H. Ge, and B. Ai., “An Efficient Nonlinear CompandingTransform for Reducing PAPR of OFDM Signals”, IEEE Trans. Broadcast.,vol. 58, no. 4, pp.: 677-684, December 2012, and the two-component,pieceweise linear compander (LIN2) in Y. Wang, J.-H. Ge, L.-H. Wang, J.Li, and B. Ai., “Nonlinear Companding Transform for Reduction ofPeak-to-Average Power ratio in OFDM Systems”, IEEE Trans. Broadcast.,vol. 59, no. 2, pp.: 369-375, June 2013. The LIN compander wasstate-of-the-art as of December 2012. The LIN compander was improved tocreate the LIN2 compander. A cutoff value of A=1.6281 was chosen, asthis appears to be the maximum cutoff value obtainable with the LIN2compander. The LIN compander can attain a maximum cutoff value of aboutA=1.44, and this value was used for the LIN compander results. Thereforethe Lagrange compander solution can provide higher cutoff values thancan be obtained with the LIN and LIN2 solutions, thereby providingadditional compander design flexibility. Two Lagrange companders wereused; the optimal solution (LG) with a perturbed a₁>0 to generate anon-negative solution, and the epsilon-perturbed compander (EP), whichtakes the optimal Lagrange compander and perturbs the tail using optimalperturbation values. For each plot, the curve labels have the followingmeaning: LG=optimal Lagrange, EP=epsilon-perturbed Lagrange, LW=LIN2,WG=LIN.

FIG. 16 contains symbol error plots for the four companders. FIG. 16shows that the LG and EP companders provide an improvement in symbolerror rate performance over the LIN and LIN2 companders.

This improvement in symbol error rate comes at the cost of a smallreduction in PAPR performance relative to the LIN2 approach, shown inFIG. 17 which contains the CCDF results. The LIN approach, under thepresent scenario, gives the best performance due to the significantlyreduced cutoff value, A=1.44 versus A=1.6281.

The Lagrange companders, however, provide a small performanceimprovement in out-of-band power rejection over the LIN2 and asignificant performance advantage over the LIN approach (FIG. 18). TheLagrange approach shows an improvement of about 0.4 dB over the LIN2approach and over 1.0 dB in out-of-band power rejection over the LINapproach.

Sample Lagrange companders, for different cutoff values, are shown inFIG. 19.

3. CONCLUSION

Those skilled in the art will appreciate that we have developed a familyof companders for PAPR reduction in OFDM signals using a constrainedoptimization approach. We derived companders and decompanders anddemonstrated that the new companders can provide performanceimprovements over current state-of-the art solutions. The new companderscan also provide solutions over ranges of the cutoff values where thecurrent state-of-the-art companders fail to exist. The set of compandersdeveloped in this paper increase the solution set of companders totradeoff between demodulation performance, PAPR reduction, andout-of-band power rejection.

The foregoing description of the embodiments of the invention has beenpresented for the purposes of illustration and description. It is notintended to be exhaustive or to limit the invention to the precise formdisclosed. Many modifications and variations are possible in light ofthis disclosure. It is intended that the scope of the invention belimited not by this detailed description, but rather by the claimsappended hereto.

What is claimed is:
 1. A signal transmitting apparatus for an orthogonalfrequency division multiplexing (OFDM) system, comprising: a companderconfigured to compress and expand a transmitted signal wherein thecompander comprises a Rayleigh amplitude distribution to a constrainedoptimization problem which may be solved using Lagrange multipliers. 2.A signal transmitting apparatus for orthogonal frequency divisionmultiplexing (OFDM) system, comprising: a companding feedback moduleconfigured to compress and expand a transmitted signal; wherein thecompanding feedback module comprises a Rayleigh amplitude distributionto a constrained optimization problem which may be solved using Lagrangemultipliers.
 3. A method of deriving a compander for orthogonalfrequency division multiplexing comprising: converting a Rayleighamplitude distribution to a constrained optimization problem which maybe solved using Lagrange multipliers.
 4. The method of claim 3 whereinsaid Rayleigh amplitude distribution is decomposed into disjoint regionsover which a plurality of amplitude weighting functions are used.
 5. Themethod of claim 3 wherein said constrained optimization approachcomprises optimizing an objective function with respect to changeableelements in the presence of constraints on those elements.
 6. The methodof claim 5 wherein said Rayleigh amplitude distribution is a probabilitydistribution function and a first said constraint is that said Rayleighamplitude distribution must integrate to unity.
 7. The method of claim 5wherein a second said constraint is that power must be constant acrossthe distribution.
 8. The method of claim 5 wherein an optimal solutionis minimally perturbed to generate only non-negative solutions.
 9. Themethod of claim 5 wherein said Rayleigh amplitude distribution isdefined as:${f_{R\;}(x)} = {\frac{2x}{\sigma^{2}}{^{- \frac{x^{2}}{\sigma^{2}}}.}}$10. The method of claim 9 wherein further comprising determining afunction g(x) which minimizes: ∫₀^(A)(f_(R)(x) − g(x))² x.
 11. Themethod of claim 10 wherein g(x) is chosen to be of piecewise linearform.
 12. The method of claim 11 wherein g(x) is defined as${g(x)} = {\sum\limits_{i = 1}^{N}{{U_{i}(x)}.}}$
 13. The method ofclaim 3 wherein the dynamics of the system are defined as: Λ(β₁, . . .β_(N+1), λ₁, λ₂)=F(β₁, . . . β_(N+1))+λ₁[g₁(β₁, . . .β_(N+1))−1]+λ₂[g₂(β₁, . . . β_(N+1))−σ²].
 14. The method of claim 13wherein a first constraint, that the cumulative distribution functionmust integrate to unity, is expressed as: g₁(β₁, . . . β_(N+1))=1. 15.The method of claim 13 wherein a second constraint, that power must beconstant, is expressed as: g₂(β₁, . . . β_(N+1))=σ².
 16. The method ofclaim 3 wherein said constrained optimization problem can be defined as:Minimize: F(β₁, …  B_(N + 1)) = ∫₀^(A) (f_(R)(x) − g(x))²x  forβ_(i), subject to ∫₀^(A) g(x)x = 1  and∫₀^(A) x²g(x)x = ∫₀^(∞)x²f_(R)(x) x = σ²,  using Lagrangemultipliers.
 17. The method of claim 3 wherein said compander is definedas:${F_{y}^{- 1}{F_{R}(x)}} = {\sqrt[ \pm ]{\frac{1 - ^{- \frac{x^{2}}{\sigma^{2}}} - Z_{i - 1}}{\gamma} - \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)} - \frac{\delta}{2\gamma}}$for a_(i)≦x<a_(i+1), wherein$Z_{i} = {Z_{i - 1} + {\frac{{\beta_{i + 1}\left( {a_{i + 1} - a_{i}} \right)} + {\beta_{i}\left( {a_{i + 1} - a_{i}} \right)}}{2}.}}$18. The method of claim 17 wherein a decompander for use with saidcompander is defined as:${x = {{{sgn}(x)}\sqrt{{- \sigma^{2}}{\ln \left( {{- {\gamma \left( {y + \frac{\delta}{2\gamma}} \right)}^{2}} - e + \frac{\delta^{2}}{4\gamma} + 1 - Z_{i - 1}} \right)}}}},$where the sgn function refers to retaining either the positive ornegative square-root, depending on the sign of x; and wherein the rangeof said compander can be defined as, when a_(i)≦x<a_(i+1), they_(MIN)≦y<y_(MAX) where: y_(MIN)=min{y(a_(i)),y(a_(i+1))},y_(MAX)=max{y(a_(i)),y(a_(i+1))}, subject to${y\left( a_{i} \right)} = {\sqrt[ \pm ]{\frac{1 - ^{- \frac{a_{i}^{2}}{\sigma^{2}}} - Z_{i - 1}}{\gamma} - \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)} - \frac{\delta}{2\gamma}}$ and${y\left( a_{i + 1} \right)} = {\sqrt[ \pm ]{\frac{1 - ^{- \frac{a_{i + 1}^{2}}{\sigma^{2}}} - Z_{i - 1}}{\gamma} - \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)} - {\frac{\delta}{2\gamma}.}}$19. The method of claim 18 wherein a Hessian Matrix is used to determinecritical point types.